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Multidimensional Unstructured-Grid Liquid Rocket Engine 
Nozzle Performance and Heat Transfer Analysis 


Ten-See Wang* 

NASA MarshallSpace Flight Center,Huntsville f Alabama t 35812 


The objectives of this study are to conduct a unified computational analysis for 
computing the design parameters such as axial thrust, convective and radiative wall heat 
fluxes for liquid rocket engine nozzles, so as to develop a computational strategy for 
computing those design parameters through parametric investigations. The computational 
methodology is based on a multidimensional, finite-volume, turbulent, chemically reacting, 
radiating, unstructured-grid, and pressure-based formulation, with grid refinement 
capabilities. Systematic parametric studies on effects of wall boundary conditions, 
combustion chemistry, radiation coupling, computational cell shape, and grid refinement 
were performed and assessed. Comparisons of the computed axial thrust performance, flow 
features, and wall heat fluxes with those of available test data and design calculations are 
presented. 


Nomenclature 

C 2f C 5 ,C>= turbulence modeling constants, 1.15, L9, 0.25, and 0.09. 

= heat capacity 
= diffusivity 
= total enthalpy 
= static enthalpy 
= radiative intensity 
= thermal conductivity 
= turbulent kinetic energy 
= nondimensional pressure 
= nondimensional heat flux 
= recovery factor 
= location coordinate 
= nondimensional temperature 
= law-of-the-wall temperature 
= time, s 

v,w = mean velocities in three directions 
= wall friction velocity 
= nondimensional Cartesian coordinates 
= turbulent kinetic energy dissipation rate 
= energy dissipation contribution 
= absorption coefficient 
= viscosity 

= turbulent eddy viscosity (=pC H k 2 /e) 

= turbulent kinetic energy production 
= density 

= turbulence modeling constants 
r = shear stress 

Q = direction vector. Q' denotes the leaving radiative intensity direction 

co = chemical species production rate 


* Staff Consultant, Applied Fluid Dynamics Analysis Group, TD64, Senior Member AIAA. 
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Subscripts 

b ~ black body 

c = convective 

cl = centerline 

p = off-wall (wall-function) point 

r = radiative 

t = turbulent flow 


w = wall 


I. Introduction 

T he two major factors in rocket engine design, performance and integrity (convective heat transfer), are often 
conducted separately. As a result, the final design based on performance may have to be altered due to 
convective heating considerations, and vice versa, resulting in delays and compromises. Recently, radiative heating 
has been generating concerns due to renewed interest in hydrocarbon engines. Those combined motivated this study 
to perform and demonstrate a unified analysis for the computation of those design parameters in a simultaneous 
fashion. Systematic parametric studies on effects of wall boundary conditions, combustion chemistry, radiation 
coupling, computational cell shape, and grid refinement were performed and assessed, in order to come up with a 
strategy for efficient and realistic computations. 

A nozzle axial force analysis, 1 and a conjugate convective heat transfer analysis 2 for axisymmetric Space Shuttle 
Main Engine (SSME) thruster were reported in the 1990’s, using a structured-grid, multi-zone computational fluid 
dynamics (CFD) solver FDNS. As the requirements for parallel computing efficiency and faster grid generation 
arise, an unstructured-grid CFD methodology UNIC was developed recently through several activities, namely the 
launch vehicle base-heating, 3 Laser propulsion, 4 and stage- separation. 5 This unstructured-grid CFD methodology is 
refined in this study to conduct the unified axial force, convective and radiative heat transfer analyses on the SSME 
thruster flowfield, simulating the hot-firing operation at sea level. Both axisymmetric and three-dimensional 
analyses were performed and a final computational strategy reported. 


II. Computational Methodology 

The time- varying transport equations of continuity, species continuity, momentum, global energy (total 
enthalpy), turbulent kinetic energy, and turbulent kinetic energy dissipation can be written as: 
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where the energy dissipation contribution 0 can be expressed as: 


0 =-. 


dxi 




and the shear stress ^ can be expressed as: 
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A predictor and corrector solution algorithm was employed to provide coupling of the fluid governing equations. 
A second-order central-difference scheme was employed to discretize the diffusion fluxes and source terms of the 
governing equations. For the convective terms, a second-order total variation diminishing difference scheme was 
used in this effort. To enhance the temporal accuracy, a second-order backward difference scheme was employed to 
discretize the temporal terms. A point-implicit (operator splitting) method was employed to solve the chemistry 
system. 

An extended k-e turbulence model 6 was used to describe the turbulence. A 7-species, 9-reaction detailed 
mechanism 7 was used to describe the finite-rate, hydrogen/oxygen (H 2 /0 2 ) afterburning chemical kinetics. The 
seven species are H 2 , 0 2 , H 2 0, O, H, OH, and N 2 . 

A modified wall function approach was employed to provide wall boundary layer solutions that are less sensitive 
to the near- wall grid spacing. Consequently, the model has combined the advantages of both the integrated-to-the- 
wall approach and the conventional law-of-the-wall approach by incorporating a complete velocity profile and a 
universal temperature profile 7 . This approach is especially useful in three-dimensional (3-D) applications. 

The convective heat transfer follows the modified Newtonian law 

Q c =(pu T /T + )[h w -h p -R(u p 2 /2} 

The radiative heat transfer is analyzed by solving the radiative transfer equation 

(Q • V)/(r, Q) = ~*7(r, Q) + Kl b (r) 

Discrete ordinate method was used to solve the radiative transfer equation and H 2 0 is the major radiating medium. 
The spectral-line based weighted-sum-of-gray-gases model 8 was used to calculate the total emissivity and 
absorptivity of the radiating medium. Details of the numerical algorithm can be found in Refs 3-5. The radiative 
heat flux is given by the integration of the wall leaving radiative intensities 


Q r = j I(r,Q )| n.Q \d£l 
n.QT < 0 


III. Computational Grid Generation 

The flowfields of four axisymmetric and four 3-D grids were computed during the course of the study. The 
results from two representative axisymmetric and two 3-D grids are reported for conciseness. These grids are hybrid 
grids and can be classified into two groups: the structure-grid dominated and the unstructured-grid dominated. 
Figure 1 shows the layout of an unstructured-grid dominated hybrid grid axl, which is an axisymmetric grid and has 
4 layers of structured (quadrilateral) grid surrounding the solid walls, while the rest of the domain filled with 
unstructured (triangular) cells. These structured- grid layers are used to ensure proper wall boundary layer 
development. The layout of a structure-cell dominated hybrid grid ax6 is shown in Fig. 2. Axisymmetric grid ax6 
has the interior of the thruster and plume region filled with quadrilateral cells, while the rest of the domain filled 
with triangular cells. The structured-grid layers used in grid axl are also embedded in grid ax6, and in 3-D grids 
3d6 and 3d9 (Figs 3-4), such that the boundary layer development for all grids is similar. Figure 3 shows the layout 
of the hybrid 3-D grid 3d6. It was constructed by rotating grid ax6 72 times for 360 degrees. Figure 4 shows the 
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layout of the hybrid grid 3d9. These computational grids were generated using the software package GRIDGEN 9 
Table 1 shows the total number of points and cells in these four grids. The structured cells in grids 3d6 and 3d9 are 
hexahedral elements. The unstructured cells in grid 3d9 are tetrahedral elements, while the unstructured cells in grid 
3d6 are prismatic elements. Note that all 4 grids were computed as a single zone, thus avoiding the interface 
complexities commonly seen in multi-zonal grids. 


Table 1. Number of points and cells of the axisymmetric and 3-D gri ds. 


Grid 

# points 

# cells 

# structured 
cells 

# unstructured 
cells 

Axl 

17,509 

30,578 

2,016 

28,562 

Ax6 

17,391 

17,710 

15,300 

2,410 

3d6 

1,286,934 

1,275,120 

1,101,600 

173,520 

3d9 

418,165 

1,732,081 

227,984 

1,504,097 


IV. Boundary Conditions and Run Matrix 

Fixed total condition was used for the inlet of the thruster and the outer boundary. A total pressure of 1 atm was 
specified for the outer boundary in order to simulate the nozzle hot-firing at sea level. No-slip boundary condition 
was specified for the thruster walls. Symmetry condition was applied to the centerline for axisymmetric cases. The 
CEC program 10 was used to obtain the chamber equilibrium species composition for use at the thruster inlet. 

The run matrix is shown in Table 2. These cases were built up systematically in order to understand the grid 
effects such as cell shape and grid refinement, and the modeling effects such as wall boundary condition, chemistry, 
and radiation. For the convenience of presentation, abbreviated letters are used to represent different cases in the 
run matrix. For example, case fz represents parametric conditions of adiabatic wall and frozen chemistry, while case 
frcgr uses parametric conditions of cooled wall and finite-rate chemistry, with radiation coupling and grid 
refinement turned on. Due to the limitation of current resources, grid refinement was not preformed for the 3-D 
cases. 


Table 2. Run matrix 


Case 

wall 

chemistry 

Grid 

refinement 

Radiation 

fz 

adiabatic 

frozen 

no 

no 

eq 

adiabatic 

equilibrium 

no 

no 

& 

adiabatic 

Finite-rate 

no 

no 

frc 

cooled 

Finite-rate 

no 

no 

frcr 

cooled 

Finite-rate 

no 

yes 

frcg 

cooled 

Finite-rate 

yes 

no 

frcgr 

cooled 

Finite-rate 

yes 

yes 


V. Results and Discussion 

The computation was performed on a cluster machine using 4 processors for each axisymmetric case and 32 
processors for each 3-D case. A global time step of 1 ps was used. Figure 5 shows a comparison of the computed 
adiabatic wall temperatures for grid ax6. Similar, if not identical wall temperature profiles were obtained for grid 
axl and are not shown. The computed wall temperature for the frozen chemistry case is constant, indicating the 
conservation laws are satisfied. Those for the equilibrium and finite-rate chemistry cases increase first after the 
throat, due to the recombination of chemical species to become H 2 0. The temperatures for those two cases then 
decrease as H 2 0 dissociates. Of interest is the temperature for the equilibrium case, it drops continuously until it 
closes to that of the frozen flow, near the nozzle exit. This is expected since the stagnation temperature is very close 
to the chamber temperature with which the frozen composition was determined with an equilibrium solution. This 
also implies the equilibrium chemistry probably dissociates the H 2 0 at too fast a rate inside the nozzle. A specified 
cooled wall temperature profile 1 is also shown in Fig. 5 which was determined through a separate conjugate heat 
transfer calculation; this temperature profile is used later as a boundary condition to consider the effect of heat loss 
to the regenerative coolant channels. 
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Figure 6 shows the computed Mach number contours for cases fz and frcgr for grids axl and ax6, respectively; 
those for other cases are similar to those of case frcgr and are not shown. These figures show the captured nozzle 
flow features (nozzle shock, lip shock, triple point, Mach disc, shock reflection, and shear layer/shock interaction). 
In general, all cases capture the flow features reasonably well, except the frozen flow case in which a curved Mach 
disk was obtained. It can also be seen that the nozzle shocks appear to be sharper in the contours of the structured- 
cell dominated grid ax6 than those of the unstructured-element dominated grid axl. The sharpest nozzle flow 
features are captured with grid refinement on grid ax6, while the added radiation changes the flow features only 
slightly. 

The significance of a curved disk is that a large flow recirculation appears behind the curved disk. The 
occurrence of the curved disk may be attributed to the difference in thermodynamics between the frozen flow and 
chemically reacting flows. As shown in the centerline H 2 0 mass fraction, specific heat ratio and Mach number 
profiles in Fig. 7, the inability of recombining the species of the frozen flow results in much higher specific heat 
ratios than those of the reacting flows, which in turn produces high Mach numbers along the centerline. The higher 
shock strength leads to higher total pressure loss across the shock, which causes the shock center to retreat and 
consequently an overall curved disk. On the other hand, the curves of equilibrium chemistry closely follow those of 
finite-rate chemistry. This is because the centerline temperatures drop continuously and are much lower than the 
chamber temperature (Fig. 8), hence the dissociation process occurring on the adiabatic wall is frozen on the 
centerline. Furthermore, the curved disk phenomenon happens both in grid axl and ax6, hence it is cell-shape 
independent and thermodynamics induced. 

Figure 8 shows a comparison of the thruster centerline temperatures for grid ax6. The frozen chemistry gives the 
lowest bound while all other cases group together as an upper bound, while the result from Ref. 1 falls in between. 
As discussed above, the low frozen chemistry curve is caused by the thermodynamics. As for the result from Ref. 1, 
it is speculated that an older thermodynamics database was used then. A comparison of the thruster wall pressures is 
shown in Fig. 9. The computed results from all cases appear to group together and agree reasonably well with the 
test data. Figure 10 shows a comparison of thruster centerline pressures. All predictions agree reasonably well, 
except for the frozen flow case that deviates lower near the nozzle lip. 

Figure 1 1 shows the computed convective heat fluxes for grid ax6. As expected, the peak convective heat fluxes 
occur at the throat (x = 0) for all cases. The refined grid gives a slightly lower peak heat flux. Radiation does not 
affect the convective heat flux, because the maximum radiative heat flux is about two orders-of-magnitude lower 
than that of convection (Fig. 12). All predictions compare reasonably well with those of the three design 
methods. 2,11 Result from Fig. 11 demonstrates that the wall boundary layers were captured reasonably well. The 
difference in the initial heat fluxes is caused by the difference in ways of initiating the boundary layers among 
different methods and the significance of which is negligible in comparison to the peak heat flux. 

Figure 12 shows the computed radiative heat flux for grid ax6 while cooled wall, finite-rate chemistry, and grid 
refinement were used as operating conditions. As expected, high radiative heating occurs inside the combustion 
chamber within which the high temperature and high H 2 0 concentration are prevalent. As the propulsive flow 
expands past the throat, the temperature drops, hence the low radiative heat flux. The peak radiative heat flux is 
about two orders of magnitude lower than that of the convective heat flux, which is reasonable for a hydrogen fueled 
engine. In current methodology, the injector faceplate is modeled as a black body. In order to compare the 
predicted radiation with that of a plume radiation code GASRAD, 12 which does not model the injector faceplate, 
another run was performed by setting the temperature of the injector faceplate to 300 deg. K, effectively turning off 
the black body radiation. The structured-grid solution from Ref. 1 was used as input for GASRAD radiation 
calculation, since GASRAD can not read unstructured-grid information. It can be seen that the result from turning 
off the blackbody radiation at the inlet, using a weighted-sum-of-gray-gases absorption model, compares reasonably 
well with that of GASRAD in which a narrow band absorption model was used. It should be pointed out that 
GASRAD reads in flow solution for a decoupled radiation solution, and current methodology solves the flow 
equations and radiative transfer equation simultaneously. The computed peak value is higher when the black body 
radiation is included at the inlet, as expected. It is also noted that GASRAD was developed for the prediction of 
plume radiation, hence it does not consider the re-radiation from the solid walls. In addition, it solves the line-of- 
sight equation and not the radiative transport equation. 


Table 4. Comparison of SSMDE thrust ch amber specific impulses (ISP’s) for axisymmetric cases 



Axl 

Ax6 

fz 

438.7 

439.7 

eq 

455.6 

456.0 

fr 

455.2 

455.6 
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frc 

452.4 

452.6 

frcg 

452.9 

453.7 

frcgr 

452.5 

453.3 

Data 

453.3 


Table 4 shows the comparison of computed SSME thrust chamber specific impulses, or the axial thrust 
performances for the axisymmetric cases. The frozen flow calculations give too low an axial force, even for the 
adiabatic wall assumption that assumes zero wall heat loss. This is again caused by inadequate heat capacity 
distributions forced by a fixed species composition. The reacting flow (with adiabatic wall) cases overpredict the 
data for about 2~3 s, with the equilibrium case giving the highest values. When the wall heat loss is considered, the 
axial force predictions become very close to the data. The quadrilateral cell dominated grid ax6 appears to predict 
slightly better ISP’s than those of the triangular cell dominated grid ax I. Within grid ax6, the grid refinement and 
radiation options offer the best agreement. 

Figure 13 shows the computed temperature contours for grid 3d6, case frc. Similar to the Mach number 
contours, the temperature contours also show the captured nozzle flow physics such as the nozzle shock, lip shock, 
triple point, Mach disc, shock reflection, and shear layer/shock interaction. Two perpendicular planes are used to 
give the Mach disc a three-dimensional feel. The high temperature in the mixing layer indicates after burning. 

Figure 14 shows a comparison of computed thruster centerline temperatures. The centerline temperature of grid 
3d6 matches that of grid ax6 reasonable well, except inside the chamber where the temperature of grid 3d6 is 
slightly lower. Figure 15 compares the wall pressures. The wall pressures of grid 3d6 and ax6 overlap and both 
compare reasonably well with the data. Figure 16 compares the centerline pressures. The centerline pressure of grid 
3d6 coincides with that of grid ax6, until the nozzle lip where the pressure of grid ax6 is slightly higher. 

Figure 17 shows a comparison of convective wall heat fluxes. The computed heat fluxes agree reasonably well 
with those of the design methods. The predictions for grid 3d6 overlap with those of grid ax6. The radiation does 
not affect the convective heat fluxes of grid 3d6, again due to the relative low radiative heat fluxes inside a H2/O2 
engine. Figure 18 shows a comparison of the computed radiative wall heat fluxes. Similar to the result of the 
axisymmetric cases (Fig. 12), the computed radiative fluxes using a weighted-sum-of-gray gases absorption model 
compares reasonably well with that of GASRAD using a narrow band absorption model, while turning off the 
blackbody radiation at the inlet. And the predicted radiative heat flux is higher while turning on the blackbody 
radiation at the inlet. 

Table 5. Comparison of SSME thrust ch amber specific impulses (ISP’s) for 3-D cases 



3d6 

3d9 

fz 

439.6 

436.5 



456.4 

453.3 

fr 

454.9 

452.5 

frc 

453.2 

450.0 

frcr 

453.1 

449.9 

Data 

453.3 


Table 5 shows the comparison of computed specific impulses for the 3-D cases. The qualitative trend among the 
cases is very similar to the corresponding axisymmetric cases. The results from grid 3d9 are consistently lower than 
those of grid 3d6. This is because the effective cell density of grid 3d9 is less than that of 3d6, although the total 
number of cells in grid 3d9 is higher than that of grid 3d6 (Table 1). As a general rule of thumb, the accuracy of two 
tetrahedral cells is approximately equivalent to that of one hexagonal cell. On the other hand, since the number of 
cells in grid 3d9 is more than those in grid 3d6, it costs more to run grid 3d9. This demonstrates that the structured- 
cell dominated grid 3d6 is favorable both in terms of accuracy and computational efficiency, similar to the findings 
in the axisymmetric cases. This also agrees with the result of Huynh’s Fourier analysis 13 that the upwind scheme 
prefers structured meshes. Within grid 3d6, again the result of case frc compares very well with that of the 
measurement, while the addition of radiation (case frcr) changes the value only slightly. 


VI. Conclusions 

Unified computational analyses for computing the design parameters such as the axial thrust, convective and 
radiative wall heat fluxes for liquid rocket engine thrusters were conducted, in order to develop a computational 
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strategy for computing those design parameters through parametric investigations. The computational methodology 
is based on a multi-dimensional, finite-volume, turbulent, chemically reacting, radiating, unstructured-grid, and 
pressure-based formulation. Systematic parametric studies on effects of wall boundary conditions, combustion 
chemistry, radiation coupling, computational cell shape, and grid refinement were performed and assessed. Under 
the computational framework of this study, it is found that the structured-mesh performed favorably than the 
unstructured-mesh. The effects of radiation coupling and grid refinement were demonstrated. Finite-rate chemistry 
performed favorably than the equilibrium chemistry, while the frozen chemistry is undesirable, due to 
thermodynamics considerations. For regeneratively cooled engines, incorporating the effect of heat loss drastically 
improves the axial force predictions. The computed flow physics, axial thrust performance, and wall heat fluxes 
compared well with those of available test data and design calculations, when the favored computational strategy 
was used. 
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Fig. 1 The layout of hybrid grid axl. Top: the overall 
grid. Bottom left: close-up near the throat Bottom 
right: close-up near the nozzle lip. 


Fig. 3 Layout of hybrid grid 3d6. Upper figure: an 
overall view. Lower left: a cross-sectional cut through 
the axis. Lower right: the exit plane. 



Fig. 2 The layout of hybrid grid ax6. Top: the overall 
grid. Bottom left: close-up near the throat Bottom 
right: close-up near the nozzle lip. 



Fig. 4 Layout of hybrid grid 3d9. Lower left insert: 
the cross-sectional cut of the thruster inlet Lower 
right insert: the cross-sectional cut of the exit plane. 
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Fig. 5 A comparison of computed and specified 
regeneratively cooled wall temperatures for grid ax6. 


Fig. 7 A comparison of computed centerline H 2 0 mass 
fractions, specific heats and Mach numbers for grid 
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Ref. 1, 1993 
Grid ax6 Casefz 
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Fig. 6 Computed Mach number contours, a) case fz, 
grid axl; b) case frcgr, grid axl; c) case fz, grid ax6; 
and d) case frcgr, grid ax6. 


Fig. 8 A comparison of thruster centerline 
temperatures for grid ax6. 
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Fig. 13 Computed temperature contours for grid 3d6, 
casefrc. 



Fig. 14 A comparison of thruster centerline 
temperatures. 


American Institute 


i 


2 


X 


3 


Fig. 15 A comparison of thruster wall pressures. 



Fig. 16 A comparison of thruster centerline pressures. 
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Fig. 17 A comparison of wall heat fluxes. Fig. 18 A comparison of radiative wall heat fluxes. 
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